---
title: "Random effects resubmission"
format: pdf
editor: visual
---

#This code lays out the necessary libraries and reads in the data.

```{r}
# Load necessary libraries
library(tidyverse)
library(readxl)
library(metafor)
library(stargazer)

data <- read_excel("C:/Users/Downloads/Black+Palestine+Solidarity_May+17,+2024_19.19.csv.xlsx")
```

#This code cleans the data and creates indexes of the solidarity variables.

```{r}
# Code demographics
data <- data %>%
  mutate(
    ages1 = as.numeric(factor(Age)),
    age = ages1 + 17,
    gender = as.numeric(factor(Gender)),
    female = ifelse(gender == 2, 1, ifelse(gender == 1, 0, NA)),
    educ = as.numeric(factor(Education)),
    partyID = as.numeric(factor(Partisanship)),
    leaners = as.numeric(factor(Ind_Lean))
  )

# Recode values
data$ages1[data$ages1 %in% c(53, 54, 55, 56, 57)] <- NA
data$educ[data$educ %in% c(5, 6)] <- NA
data$partyID[data$partyID %in% c(6, 7)] <- NA
data$leaners[data$leaners %in% c(4, 5)] <- NA

# Create manipulation variable
library(dplyr)

data <- data %>%
  mutate(
    control = as.numeric(factor(Control_Turtle)),
    SharedDisc1 = as.numeric(factor(SharedDisc)),
    shareddisc = ifelse(SharedDisc1 == 1, 1, ifelse(control == 1, 0, NA))
  )

data <- data %>%
  mutate(
    shareddisc = case_when(
      Control_Turtle == 1 ~ 0,  # Assign 0 for control group
      SharedDisc == 1 ~ 1,  # Assign 1 for treatment group
      TRUE ~ NA_real_  # Assign NA if neither condition is true
    )
  )

data <- data %>% filter(!is.na(shareddisc))
# Code solidarity variables
data <- data %>%
  mutate(
    bondP = as.numeric(factor(Bond_Palestine)),
    alliesP = as.numeric(factor(Allies_Palestine)),
    fateP = as.numeric(factor(CFate_Palestine))
  )
data$bondP[data$bondP %in% c(6, 7)] <- NA
data$alliesP[data$alliesP %in% c(6, 7)] <- NA
data$fateP[data$fateP %in% c(6, 7)] <- NA

# Create solidarity indexes
data <- data %>%
  mutate(
    solidarity1 = (bondP + alliesP + fateP) / 3,
    solidarity = (solidarity1 - 1) / 4
  )

# Regression models
summary(lm(solidarity ~ shareddisc, data = data))

# Code policy items
data <- data %>%
  mutate(
    policy1 = as.numeric(factor(Policy1)),
    policy2 = as.numeric(factor(Policy2)),
    policy3 = as.numeric(factor(Policy3)),
    policy4 = as.numeric(factor(Policy4))
  )
data$policy1[data$policy1 %in% c(6, 7)] <- NA
data$policy2[data$policy2 %in% c(6, 7)] <- NA
data$policy3[data$policy3 %in% c(6, 7)] <- NA
data$policy4[data$policy4 %in% c(6, 7)] <- NA

# Reverse scale for policy items
data$policy1 <- recode(data$policy1, `1` = 5, `2` = 4, `3` = 3, `4` = 2, `5` = 1)
data$policy2 <- recode(data$policy2, `1` = 5, `2` = 4, `3` = 3, `4` = 2, `5` = 1)
data$policy4 <- recode(data$policy4, `1` = 5, `2` = 4, `3` = 3, `4` = 2, `5` = 1)

# Create pro-Palestinian policy index
data <- data %>%
  mutate(
    propal1 = (policy1 + policy2 + policy3 + policy4) / 4,
    propal = (propal1 - 1) / 4
  )


```

#The code below gives the mediation pathways of the treatment on outcome, treatment on mediator, and mediator on outcome.

```{r}
data <- na.omit(data[, c("propal", "solidarity", "shareddisc")])

# Model 1: Treatment to Outcome (propal ~ shareddisc)
model_t_o <- lm(propal ~ shareddisc, data = data)
summary(model_t_o)
# Model 2: Treatment to Mediator (solidarity ~ shareddisc)
model_t_m <- lm(solidarity ~ shareddisc, data = data)
summary(model_t_m)
# Model 3: Mediator to Outcome (propal ~ solidarity)
model_m_o <- lm(propal ~ solidarity, data = data)
summary(model_m_o)

# Extract coefficients from models
beta_t_o <- summary(model_t_o)$coefficients["shareddisc", "Estimate"]
beta_t_m <- summary(model_t_m)$coefficients["shareddisc", "Estimate"]
beta_m_o <- summary(model_m_o)$coefficients["solidarity", "Estimate"]

p_value_t_o <-summary(model_t_o)$coefficients["shareddisc", "Pr(>|t|)"]
p_value_t_m <- summary(model_t_m)$coefficients["shareddisc", "Pr(>|t|)"]
p_value_m_o <- summary(model_m_o)$coefficients["solidarity", "Pr(>|t|)"]

# Print coefficients
cat("Beta for T to O: ", beta_t_o, "\n")
cat("Beta for T to M: ", beta_t_m, "\n")
cat("Beta for M to O: ", beta_m_o, "\n")

#Get the standard deviation
sd_propal <- sd(data$propal)
sd_solidarity <- sd(data$solidarity)

# Compute Cohen's d values
cohen_d_t_o <- beta_t_o / sd_propal
cohen_d_t_m <- beta_t_m / sd_solidarity
cohen_d_m_o <- beta_m_o / sd_propal

# Print Cohen's d values
cat("Cohen's d for T to O: ", cohen_d_t_o, "\n")
cat("Cohen's d for T to M: ", cohen_d_t_m, "\n")
cat("Cohen's d for M to O: ", cohen_d_m_o, "\n")

# Add Study 6 values to the study_data data frame
study_6 <- data.frame(
  study = "Study 6",
  d = c(cohen_d_t_o, cohen_d_t_m, cohen_d_m_o),  # Cohen's d values
  p_value = c(p_value_t_o, p_value_t_m, p_value_m_o),  # p-values
  n = rep(852, 3)  # Sample size for Study 6
)
#All values below taken from Table 2 of "Taking Stock of Solidarity between People of Color: A Mini Meta-Analysis of Five Experiments"

study_5 <- data.frame(
  study = "Study 5",
  d = c(0.209),  # Cohen's d values
  p_value = c(0.003),  # p-values
  n = rep(818)  # Sample size 
)
study_4 <- data.frame(
  study = "Study 4",
  d = c(0.106),  # Cohen's d values
  p_value = c(0.137),  # p-values
  n = rep(802)  # Sample size 
)
study_3 <- data.frame(
  study = "Study 3",
  d = c(0.362),  # Cohen's d values
  p_value = c(0.000),  # p-values
  n = rep(424)  # Sample size 
)
study_2 <- data.frame(
  study = "Study 2",
  d = c(0.012),  # Cohen's d values
  p_value = c(0.898),  # p-values
  n = rep(546)  # Sample size 
)
study_1 <- data.frame(
  study = "Study 1",
  d = c(0.119),  # Cohen's d values
  p_value = c(0.158),  # p-values
  n = rep(550)  # Sample size 
)

# Combine with the previous data
study_data <- rbind(study_1, study_2, study_3, study_4, study_5, study_6)

# Calculate the standard error of the effect sizes (Cohen's d) using the formula
study_data$se_d <- sqrt((1 / study_data$n) + (study_data$d^2 / (2 * (study_data$n - 1))))

# Convert Cohen's d to Fisher's Z
study_data$z <- atanh(study_data$d)

# Calculate the standard error of the Fisher's Z values
study_data$se_z <- 1 / sqrt(study_data$n - 3)

# Perform random-effects meta-analysis using the inverse-variance method on Fisher's Z
random_effects_model <- rma(yi = z, sei = se_z, data = study_data, method = "REML")

# Print the results
summary(random_effects_model)

# Convert the average Fisher's Z back to Cohen's d
weighted_avg_z <- random_effects_model$b
weighted_avg_d <- tanh(weighted_avg_z)

cat("Meta-analyzed Cohen's d: ", weighted_avg_d, "\n")
```

##Produce tables and figures

#This codes figure 3 in the paper. #Figure 3 - Effects of Shared Discrimination on PoC Solidarity, T on M

```{r}
# Define the data
study_data2 <- data.frame(
  study = c("Study 1", "Study 2", "Study 3", "Study 4", "Study 5", "Study 6"),
  d = c(0.218, 0.217, 0.242, 0.129, 0.131, -.027),  # Cohen's d values
  p_value = c(0.011, 0.013, 0.014, 0.065, 0.059, .693),  # p-values
  n = c(557, 542, 424, 804, 818, 851)  # Sample sizes
)

# Calculate the standard error of the effect sizes (Cohen's d) using the formula
study_data2$se_d <- sqrt((1 / study_data2$n) + (study_data2$d^2 / (2 * (study_data2$n - 1))))

# Convert Cohen's d to Fisher's Z
study_data2$z <- atanh(study_data2$d)

# Calculate the standard error of the Fisher's Z values
study_data2$se_z <- 1 / sqrt(study_data2$n - 3)

# Perform random-effects meta-analysis using the inverse-variance method on Fisher's Z
random_effects_model2 <- rma(yi = z, sei = se_z, data = study_data2, method = "REML")

# Print the results
summary(random_effects_model2)

# Convert the average Fisher's Z back to Cohen's d
weighted_avg_z <- random_effects_model2$b
weighted_avg_d <- tanh(weighted_avg_z)

cat("Meta-analyzed Cohen's d: ", weighted_avg_d, "\n")



forest_plot2 <- forest(random_effects_model2, 
                      slab = study_data2$study, 
                      xlab = "Cohen's d", 
                      mlab = "Meta-analyzed Cohen's d", 
                      psize = 1.5, 
                      cex = 0.8,
                      main = "Figure 3 - Effects of Shared Discrimination on PoC Solidarity")

text(x = study_data2$d, 
     y = seq(length(study_data2$d), 1, by = -1) + 0.005,  # Increase shift to fully move above
     labels = round(study_data2$d, 3), 
     pos = 3,  # Ensure labels are above the points
     cex = 0.8)

 random_effects_model2$pval  # p-value for the meta-analysis
 random_effects_model2$se        # Standard error of the meta-analysis estimate
 sum(study_data2$n) 
 
```

#Below is code for Meta results for table 3.

```{r}

#Effects of Shared Discrimination on PoC Solidarity, T on M
# Add Meta-Analyzed Effect
meta_row2 <- data.frame(
  Study = "Meta-Analyzed Effect",
  Effect_Size = round(random_effects_model2$b, 3),
  SE = round(random_effects_model2$se, 3),
  p_value = round(random_effects_model2$pval, 3)
)
print(meta_row2)

# Define the study data
study_data2 <- data.frame(
  Study = c("Study 1", "Study 2", "Study 3", "Study 4", "Study 5", "Study 6", "Meta-analyzed"),
  Effect_Size = c(0.218, 0.217, 0.242, 0.129, 0.131, -0.027, .151),  # Cohen's d values
  SE = c(0.043, 0.043, 0.05, 0.035, 0.035, 0.034, .042),  # Standard errors
  p_value = c(0.011, 0.013, 0.014, 0.065, 0.059, 0.693, .0003)  # p-values
)

study_data2$p_value[study_data2$Study == "Meta-analyzed"] <- format(0.000, digits = 3, nsmall = 3)

# Print combined data frame
print(study_data2)

# Create the table using stargazer and specify column names
stargazer(study_data2, type = "html", summary = FALSE, 
          title = "Table 3 - Effects of Shared Discrimination on PoC Solidarity", 
          out = "meta_analysis_table_3.html")

# Extract meta-analysis diagnostics
tau2 <- random_effects_model2$tau2  # Between-study variance
i2 <- random_effects_model2$I2  # Percentage of variance due to heterogeneity
q_stat <- random_effects_model2$QE  # Q-statistic for heterogeneity test
p_q <- random_effects_model2$QEp  # p-value for Q-test

# Print diagnostic statistics
cat("Meta-Analysis Heterogeneity Diagnostics:\n")
cat("Tau² (Between-study variance):", round(tau2, 4), "\n")
cat("I² (Heterogeneity %):", round(i2, 2), "%\n")
cat("Q-statistic:", round(q_stat, 3), " (p =", round(p_q, 3), ")\n")



#Effects of PoC Solidary on Pro-Outgroup Policy
meta_row_r <- data.frame(
  Study = "Meta-Analyzed (d)",  # Label for the last row
  Effect_Size = 0.778,          # Cohen's r value
  SE = NA,                     # No standard error for meta-analysis
  p_value = NA                 # No p-value for meta-analysis
)

print(meta_row_r)

meta_row_rconfirm <- data.frame(
  Study = "Meta-Analyzed Effect",
  Effect_Size = round(random_effects_model_r$b, 3),
  SE = round(random_effects_model_r$se, 3),
  p_value = round(random_effects_model_r$pval, 3)
)

print(meta_row_rconfirm)


# Define the study data
study_data_r <- data.frame(
  Study = c("Study 1", "Study 2", "Study 3", "Study 4", "Study 5", "Study 6", "Meta-analyzed (r)"),
  Effect_Size = c(0.334, 0.389, 0.414, 0.236, 0.462, 0.335, .380), # Cohen's d values
  SE = c(0.043, 0.043, 0.05, 0.035, 0.035, 0.034, .038),  # Standard errors
  p_value = c(0.001, 0.001, 0.001, 0.001, 0.001, 1.01e-35, 0.000))

study_data_r$p_value[study_data_r$Study == "Study 6"] <- format(0.000, digits = 3, nsmall = 3)
study_data_r$p_value[study_data_r$Study == "Meta-analyzed (r)"] <- format(0.000, digits = 3, nsmall = 3)

# Combine study data and meta row
study_data_combined_r <- rbind(study_data_r, meta_row_r)

# Print combined data frame
print(study_data_combined_r)

# Create the table using stargazer and specify column names
stargazer(study_data_combined_r, type = "html", summary = FALSE, 
          title = "Table 4 - Effects of PoC Solidary on Pro-Outgroup Policy", 
          out = "meta_analysis_table_r.html")

# Extract meta-analysis diagnostics
tau2 <- random_effects_model_r$tau2  # Between-study variance
i2 <- random_effects_model_r$I2  # Percentage of variance due to heterogeneity
q_stat <- random_effects_model_r$QE  # Q-statistic for heterogeneity test
p_q <- random_effects_model_r$QEp  # p-value for Q-test

# Print diagnostic statistics
cat("Meta-Analysis Heterogeneity Diagnostics:\n")
cat("Tau² (Between-study variance):", round(tau2, 4), "\n")
cat("I² (Heterogeneity %):", round(i2, 2), "%\n")
cat("Q-statistic:", round(q_stat, 3), " (p =", round(p_q, 3), ")\n")


```

#This code produces Figure 4 in the paper. #Figure 4 - PoC Solidary on Pro-Outgroup Policy, M on O

```{r}
# Define the data
study_data_r <- data.frame(
  study = c("Study 1", "Study 2", "Study 3", "Study 4", "Study 5", "Study 6"),
  r = c(0.334, 0.389, 0.414, 0.236, 0.462, .335),  # Correlation coefficients (r values)
  p_value = c(0.001, 0.001, 0.001, 0.001, 0.001, 1.01e-35),  # p-values
  n = c(547, 546, 417, 802, 818, 851)  # Sample sizes
)

# Convert correlation (r) to Fisher's Z
study_data_r$z <- atanh(study_data_r$r)

# Calculate the standard error of the Fisher's Z values
study_data_r$se_z <- 1 / sqrt(study_data_r$n - 3)

# Perform random-effects meta-analysis using the inverse-variance method on Fisher's Z
library(metafor)
random_effects_model_r <- rma(yi = z, sei = se_z, data = study_data_r, method = "REML")

# Print the results of the meta-analysis
summary(random_effects_model_r)
random_effects_model_r$pval

# Convert the average Fisher's Z back to correlation (r)
weighted_avg_z <- random_effects_model_r$b
weighted_avg_r <- tanh(weighted_avg_z)


cohens_d <- (2 * weighted_avg_r) / sqrt(1 - weighted_avg_r^2)

# Print the results
cat("Meta-analyzed correlation (r): ", weighted_avg_r, "\n")
cat("Meta-analyzed Cohen's d: ", cohens_d, "\n")


# Generate the forest plot for the correlation estimates
forest_plot_r <- forest(random_effects_model_r, 
                        slab = study_data_r$study, 
                        xlab = "Correlation Coefficient (r)", 
                        mlab = "Meta-analyzed Correlation (r)", 
                        psize = 1.5, 
                        cex = 0.8,
                        main = "Figure 4 - Effects of PoC Solidary on Pro-Outgroup Policy",
                        xlim = c(-.5, 1),
                        alim = c(0,1))
text(x = study_data_r$r, 
     y = seq(length(study_data_r$r), 1, by = -1) + 0.005,  # Increase shift to fully move above
     labels = round(study_data_r$r, 3), 
     pos = 3,  # Ensure labels are above the points
     cex = 0.8)


random_effects_model_r$pval  # p-value for the meta-analysis
 random_effects_model_r$se        # Standard error of the meta-analysis estimate
 sum(study_data_r$n) 
```

#Code below looks at diagnostic stats for figure 3.

```{r}

meta_row_r <- data.frame(
  Study = "Meta-Analyzed (d)",  # Label for the last row
  Effect_Size = 0.778,          # Cohen's r value
  SE = NA,                     # No standard error for meta-analysis
  p_value = NA                 # No p-value for meta-analysis
)

print(meta_row_r)

meta_row_rconfirm <- data.frame(
  Study = "Meta-Analyzed Effect",
  Effect_Size = round(random_effects_model_r$b, 3),
  SE = round(random_effects_model_r$se, 3),
  p_value = round(random_effects_model_r$pval, 3)
)

print(meta_row_rconfirm)


# Define the study data
study_data_r <- data.frame(
  Study = c("Study 1", "Study 2", "Study 3", "Study 4", "Study 5", "Study 6", "Meta-analyzed (r)"),
  Effect_Size = c(0.334, 0.389, 0.414, 0.236, 0.462, 0.335, .380), # Cohen's d values
  SE = c(0.043, 0.043, 0.05, 0.035, 0.035, 0.034, .038),  # Standard errors
  p_value = c(0.001, 0.001, 0.001, 0.001, 0.001, 1.01e-35, 0.000))

study_data_r$p_value[study_data_r$Study == "Study 6"] <- format(0.000, digits = 3, nsmall = 3)
study_data_r$p_value[study_data_r$Study == "Meta-analyzed (r)"] <- format(0.000, digits = 3, nsmall = 3)

# Combine study data and meta row
study_data_combined_r <- rbind(study_data_r, meta_row_r)

# Print combined data frame
print(study_data_combined_r)

# Create the table using stargazer and specify column names
stargazer(study_data_combined_r, type = "html", summary = FALSE, 
          title = "Table 4 - Effects of PoC Solidary on Pro-Outgroup Policy", 
          out = "meta_analysis_table_r.html")

# Extract meta-analysis diagnostics
tau2 <- random_effects_model_r$tau2  # Between-study variance
i2 <- random_effects_model_r$I2  # Percentage of variance due to heterogeneity
q_stat <- random_effects_model_r$QE  # Q-statistic for heterogeneity test
p_q <- random_effects_model_r$QEp  # p-value for Q-test

# Print diagnostic statistics
cat("Meta-Analysis Heterogeneity Diagnostics:\n")
cat("Tau² (Between-study variance):", round(tau2, 4), "\n")
cat("I² (Heterogeneity %):", round(i2, 2), "%\n")
cat("Q-statistic:", round(q_stat, 3), " (p =", round(p_q, 3), ")\n")



```

#Code below looks at diagnostic stats for Figure 4.

```{r}

# Extract meta-analysis diagnostics
tau2 <- random_effects_model$tau2  # Between-study variance
i2 <- random_effects_model$I2  # Percentage of variance due to heterogeneity
q_stat <- random_effects_model$QE  # Q-statistic for heterogeneity test
p_q <- random_effects_model$QEp  # p-value for Q-test

# Print diagnostic statistics
cat("Meta-Analysis Heterogeneity Diagnostics:\n")
cat("Tau² (Between-study variance):", round(tau2, 4), "\n")
cat("I² (Heterogeneity %):", round(i2, 2), "%\n")
cat("Q-statistic:", round(q_stat, 3), " (p =", round(p_q, 3), ")\n")

# Extract meta-analysis diagnostics
tau2 <- random_effects_model2$tau2  # Between-study variance
i2 <- random_effects_model2$I2  # Percentage of variance due to heterogeneity
q_stat <- random_effects_model2$QE  # Q-statistic for heterogeneity test
p_q <- random_effects_model2$QEp  # p-value for Q-test

# Print diagnostic statistics
cat("Meta-Analysis Heterogeneity Diagnostics:\n")
cat("Tau² (Between-study variance):", round(tau2, 4), "\n")
cat("I² (Heterogeneity %):", round(i2, 2), "%\n")
cat("Q-statistic:", round(q_stat, 3), " (p =", round(p_q, 3), ")\n")

#Table 4
# Extract meta-analysis diagnostics
tau2 <- random_effects_model_r$tau2  # Between-study variance
i2 <- random_effects_model_r$I2  # Percentage of variance due to heterogeneity
q_stat <- random_effects_model_r$QE  # Q-statistic for heterogeneity test
p_q <- random_effects_model_r$QEp  # p-value for Q-test

# Print diagnostic statistics
cat("Meta-Analysis Heterogeneity Diagnostics:\n")
cat("Tau² (Between-study variance):", round(tau2, 4), "\n")
cat("I² (Heterogeneity %):", round(i2, 2), "%\n")
cat("Q-statistic:", round(q_stat, 3), " (p =", round(p_q, 3), ")\n")

#Make into stargazer
# Extracting diagnostics from each model
# Table 2: Random effects model 1
tau2_2 <- random_effects_model$tau2  # Between-study variance
i2_2 <- random_effects_model$I2  # Percentage of variance due to heterogeneity
q_stat_2 <- random_effects_model$QE  # Q-statistic for heterogeneity test
p_q_2 <- random_effects_model$QEp  # p-value for Q-test

# Table 3: Random effects model 2
tau2_3 <- random_effects_model2$tau2  # Between-study variance
i2_3 <- random_effects_model2$I2  # Percentage of variance due to heterogeneity
q_stat_3 <- random_effects_model2$QE  # Q-statistic for heterogeneity test
p_q_3 <- random_effects_model2$QEp  # p-value for Q-test

# Table 4: Random effects model 3 (Cohen's r)
tau2_4 <- random_effects_model_r$tau2  # Between-study variance
i2_4 <- random_effects_model_r$I2  # Percentage of variance due to heterogeneity
q_stat_4 <- random_effects_model_r$QE  # Q-statistic for heterogeneity test
p_q_4 <- random_effects_model_r$QEp  # p-value for Q-test

# Creating a data frame to hold all diagnostics
meta_diagnostics <- data.frame(
  Model = c("Shared Discrimination on Pro-Outgroup Policy", "Shared Discrimination on PoC Solidarity", "PoC Solidary on Pro-Outgroup Policy"),
  Tau2 = c(round(tau2_2, 4), round(tau2_3, 4), round(tau2_4, 4)),
  I2 = c(round(i2_2, 2), round(i2_3, 2), round(i2_4, 2)),
  Q_Statistic = c(round(q_stat_2, 3), round(q_stat_3, 3), round(q_stat_4, 3)),
  p_value = c(round(p_q_2, 3), round(p_q_3, 3), round(p_q_4, 3))
)


meta_diagnostics$p_value[meta_diagnostics$Model == "Shared Discrimination on Pro-Outgroup Policy"] <- format(0.000, digits = 3, nsmall = 3)
meta_diagnostics$p_value[meta_diagnostics$Model == "Shared Discrimination on PoC Solidarity"] <- format(0.000, digits = 3, nsmall = 3)
meta_diagnostics$p_value[meta_diagnostics$Model == "PoC Solidary on Pro-Outgroup Policy"] <- format(0.000, digits = 3, nsmall = 3)

# Print the meta diagnostics data frame to verify
print(meta_diagnostics)

# Use stargazer to generate the table
stargazer(meta_diagnostics, type = "html", summary = FALSE, 
          title = "Meta-Analysis Heterogeneity Diagnostics", 
          out = "meta_analysis_diagnostics_table.html")

```
